syms a h rho r theta z M real
x = r*cos(theta); y=r*sin(theta);
pos = [x,y,z];
mass = int(int(int(rho*r,r,[0,a*(1-z/h)]),theta,[0,2*pi]),z,[0,h]);
rG = simplify(int(int(int(pos*rho*r,r,[0,a*(1-z/h)]),theta,[0,2*pi]),z,[0,h])/mass);
d = pos-rG;
integrand = (dot(d,d)*eye(3) - d'*d)*rho*r;
Ig = M*simplify(int(int(int(integrand,r,[0,a*(1-z/h)]),theta,[0,2*pi]),z,[0,h])/mass)
integrand = (dot(pos,pos)*eye(3) - pos'*pos)*rho*r;
Io = M*simplify(int(int(int(integrand,r,[0,a*(1-z/h)]),theta,[0,2*pi]),z,[0,h])/mass)